require "lfs"
require "unity"
require "MPIPond"
local MPIModel = require "MPIModel"
--local RCSPost = require "RCSPost"
local RCSPost = require"RCS.luaRCSPost"
local nl = require "luaNumeric"

-----------------------------------------------------------------------------
--初始化参数
-- RCS={
	-- num_threads = 10;
	-- dext = 8;                                    -- 物体到空气盒子
	-- farbox_info = { freq={1,2,3,4}; dfarbox = 6; };
	-- source_info = { freq_min = 5; freq_max = 10;  dsrc  = 3; type = "HH", theta = 45; };
	-- boundary = { xm=UPML; };
	-- units = { freq = GHz; dim = mm };
	-- mesh_freq = 12;
	-- name = "unname";
	-- save2file = false;
	-- camera = "3377.21 -3356.18 2213.26 0.4964 0.205616 0.322752 0.779192";
-- };
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
--后处理
--从模型读取数据后，进行单站RCS提取
function RCSPost:get_mono_RCS ()
	--print('self.theta', self.theta)
	
	--theta是俯仰角，长度是所有角度数据长度
	if self.theta == nil or #self.theta < 1 then return false end
	local n = #self.theta;
	
	--phi的索引是所有角度数据的中间数值
	local idx = math.floor((n+1)/2);
	self.scan_angle_gap = self.phi[idx];
	
	if not self.etheta_data or #self.etheta_data < 1 then
		return false;
	end
	
	--print('self.freq', self.freq)
	
	
	self.etheta = {}
	self.ephi = {}
	self.monoRCS = {} -- Ethotal
	  for k,v in pairs(self.freq) do
		-- self.etheta[k] = nl.log10( 4.0*pi * nl.abs(self.etheta_data[k][ idx] ) * nl.abs(self.etheta_data[k][ idx] ));
		self.etheta[k] = self.etheta_data[k][ idx]  ;
		self.ephi[k] = self.ephi_data[k][ idx]  ;
		-- self.monoRCS[k] = 10*nl.log10( 4.0*pi * nl.abs(self.etheta_data[k][ idx] ) * nl.abs(self.etheta_data[k][ idx] ));
		self.monoRCS[k] = 10*nl.log10( 4.0*pi * ( nl.abs(self.etheta_data[k][ idx] ) * nl.abs(self.etheta_data[k][ idx] ) +
								         nl.abs(self.ephi_data[k][ idx] ) * nl.abs(self.ephi_data[k][ idx] ))  );
		--print( self.etheta[ k], self.monoRCS[k])
	  end
	  return true;
end
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
--后处理
--读取文件数据后，进行单站RCS提取
function RCSPost:get_mono_RCS_from_file ()
	--print('self.theta', self.theta)
	
	--theta是俯仰角，长度是所有角度数据长度
	if self.theta == nil or #self.theta < 1 then return false end
	local n = #self.theta;
	
	--phi是方位角，索引是所有角度数据的中间数值
	local idx = math.floor((n+1)/2);
	self.scan_angle_gap = self.phi[idx];
	self.scan_angle = self.phi[idx];
		
	self.monoRCS = {} -- RCS thotal	
	
	for k,v in pairs(self.freq_idx) do
		-- print("freq:"..v);
		local data_idx = self:generate_idx(v,self.theta[1],self.phi[idx],"Etheta") 
		--print( self:generate_idx(v,self.theta[idx],self.phi[1],"ETheta") )
		local strEtheta = self.listdata[data_idx];
		if etheta ~= nil then
			local etheta =  tonumber( self.listdata[data_idx] );
			local et = nl.pow(10,etheta/10);
			local data_idx = self:generate_idx(v,self.theta[1],self.phi[idx],"Ephi") 
			local ephi =  tonumber( self.listdata[data_idx] );
			local ep = nl.pow(10,ephi/10);
			local etotal = et+ep;
			self.monoRCS[k] = 10*nl.log10(etotal); 
		else
			local data_idx = self:generate_idx(v,self.theta[1],self.phi[idx],"RCS-Total") 
			local rcstotal = tonumber( self.listdata[data_idx] );
			assert(rcstotal, 'The data sheet does not contain RCS results(导入文件不包含RCS计算结果)');
			self.monoRCS[k] = rcstotal; 
		end
		-- self.monoRCS[k] = 10*nl.log10(et); 
	end
	-- print_table(self.freq_idx)
	-- print_table(self.monoRCS);

	return true;
end
-----------------------------------------------------------------------------

local RCS = MPIModel:new( { 
	dext = 8;                                    -- 物体到空气盒子
	farbox_info = { freq={1,2,3,4}; dfarbox = 6; };
	source_info = { freq_min = 5; freq_max = 10;  dsrc  = 3; type = "HH", theta = 45; };
	boundary = { xm=UPML; };
	num_pml = 10;
	units = { freq = GHz; dim = mm };
	mesh_freq = 12;
	name = "unname";
	save2file = false;
	camera = "3377.21 -3356.18 2213.26 0.4964 0.205616 0.322752 0.779192";
	scan_angle = 0.0;
})

-----------------------------------------------------------------------------
--创建表和创建元表
-- function RCS:new (o)
	-- o = o or {};
	-- self.__index = self;
	-- setmetatable(o,self)
	-- return o;
-- end
-----------------------------------------------------------------------------
function RCS:init_results( )
    self.results = RCSPost:new();
end
-----------------------------------------------------------------------------
--打印函数
function RCS:print( )
	print( "swarm_size:"..self.swarm_size);
	print( "iterations:"..self.iterations);
	print( "inertia:"..self.inertia);
	print( "correction_factor:"..self.correction_factor);
	print( "max_val:"..self.max_val);
	print( "num_iter:"..self.num_iter);
	print( "swarm arr:"..#self.swarm)
	print( "init_type:"..#self.init_type)
	print( "parameter size:"..#self.par)
	print( "parameters:")
	print_table_table(self.par)
end
-----------------------------------------------------------------------------

function RCS:init()
	print("collectgarbage1");
	collectgarbage()
	collectgarbage()
	print( "Geometry clear all objects")
	Geometry.clear_all();
	
	--print( "Set path")
	--Geometry.set_path( self.path .. '/'.. self.name ..'.gap' );

	print( "Create parameters")	
	self:create_parameters();
	self:set_threads();
	self:set_boundary();
	self:set_units();
	self:set_mesh_freq();
	self:create_geometry();
	self:create_source();
	self:create_farbox();
	self:create_grabber();
end

function RCS:create_parameters()
	Geometry.Parameters.add_variable("incident_angle","0","方位角: zx平面的夹角")
	Geometry.Parameters.add_variable("polar_angle","0",'极化角,垂直极化VV[0],水平极化HH[90]')
	Geometry.Parameters.add_variable("dfarbox","10",'远场盒子与几何结构距离')
	Geometry.Parameters.add_variable("dext","20",'入射源与几何结构距离')
	Geometry.Parameters.add_variable("elevation_angle",tostring(self.elevation_angle),'俯仰角:与xy平面的夹角')	
end

function RCS:create_scaters()
	local brick = Geometry.Brick:new();
	brick:pt1(0,0,0);
	brick:pt2(10,10,10);
	
	bbox = brick:bbox()

	return {brick}, bbox
end

function RCS:create_geometry()
	print( "create geometry...")
	local dfarbox = self.dfarbox_ --or 0.1
	local dext = self.dext-- or 0.1
	
	objlist,self.bbox = self:create_scaters();	
	-- self.bbox:extend(dext,dext,dext)
	-- self.bbox.bbox_x1 = self.bbox.bbox_x1 - dext ;
	-- self.bbox.bbox_y1 = self.bbox.bbox_y1 - dext;
	-- self.bbox.bbox_z1 = self.bbox.bbox_z1 - dext;
	-- self.bbox.bbox_x2 = self.bbox.bbox_x2 + dext;
	-- self.bbox.bbox_y2 = self.bbox.bbox_y2 + dext;
	-- self.bbox.bbox_z2 = self.bbox.bbox_z2 + dext;
	
	 ---------------------------
	--Boundary Box
	print( "    create boundary box")
	local brick = Geometry.Brick:new()
	brick :name("air box")
	-- brick :pt1(self.bbox:xmin(),self.bbox:ymin(),self.bbox:zmin())
	-- brick :pt2(self.bbox:xmax(),self.bbox:ymax(),self.bbox:zmax())
	
	brick :pt1(self.bbox:xmin()..'-dext',self.bbox:ymin()..'-dext',self.bbox:zmin()..'-dext')
	brick :pt2(self.bbox:xmax()..'+dext',self.bbox:ymax()..'+dext',self.bbox:zmax()..'+dext')
	
	brick :outline_width(2)
	brick :is_fill(false)
	Geometry.AddShape(brick)
	----------------------------

         for k,v in pairs ( objlist ) do
		print( "    create object: ",v:name())
		Geometry.AddShape(v)
	end
end


function RCS:create_source(  )
	print( "create srouce")	
	-------------------------------------
	local polar = self.source_info.polarization;
	local freqmin = self.source_info.freq_min; --����Դ��ʼƵ��
	local freqmax = self.source_info.freq_max; --����Դ��ֹƵ��
	local dsrc = self.source_info.dsrc or 0;
	local theta = self.source_info.theta;
	assert(freqmin,"freqmin is not set in source function(.source_info.freqmin_)");
	assert(freqmax,"freqmax is not set in source function(.source_info.freqmax_)");
	assert(type,"src type is not set in source function(self.type_)");
	-------------------------------------

	-------------------------------------
	--Pulse
	local signal_gauss = Geometry.NewGaussPulse("Gauss Pulse")
	signal_gauss :freq_low(freqmin)
	signal_gauss :freq_high(freqmax)
	self.psig = signal_gauss
	-------------------------------------


	-------------------------------------
	--Source
	local src= Geometry.PlaneWaveSrc:new();
	if polar == "HH" then                    
	    -- src :name("HH PW src")
	    src :polarization_angle(90)
	else
	   -- src :name("VV PW src")
	end
	
	src :polarization_angle('polar_angle')

	--if polar == 'HH' then
		--Geometry.Parameters.set_variable("polar_angle","90")
	--else
		--Geometry.Parameters.set_variable("polar_angle","0")
	--end
	
	src :pt1(self.bbox:xmin()..'-dext', self.bbox:ymin()..'-dext', self.bbox:zmin()..'-dext')
	src :pt2(self.bbox:xmax()..'+dext', self.bbox:ymax()..'+dext', self.bbox:zmax()..'+dext')
	
	src :signal(self.psig)
	src :direction(forward)
	src :point_size(2)

	src:enable_plane(yx)
	src:enable_plane(yz)
	src:enable_plane(zy)
	src:enable_plane(xz)
	src:enable_plane(zx)
	--src :is_point(false)
	 -- src :is_fill(false)
	--src :is_OutLine(true)
	src:rotate(yAxis,'-elevation_angle-90')	
	src:rotate(zAxis,'incident_angle')
	src:translate( self.bbox:xmin()+self.bbox:xsize()/2,
			     self.bbox:ymin()+self.bbox:ysize()/2,
			     self.bbox:zmin()+self.bbox:zsize()/2 )
	src :outline_width(1.5)
	Geometry.AddSrc(src)
	-------------------------------------
end


function RCS:create_farbox(  )
	print( "create farbox")	
	local dfarbox = self.farbox_info.dfarbox 
	local farbox_info = self.farbox_info;
	local strFreq="";
	for k,v in pairs (self.farbox_info.freq) do
		strFreq = strFreq..v..";"
	end
	---------------------------
	--FarField Box
	local farbox = Geometry.Near2FarFieldBox:new("fbox")
	farbox :pt1(self.bbox:xmin()+dfarbox,self.bbox:ymin()+dfarbox,self.bbox:zmin()+dfarbox)
	farbox :pt2(self.bbox:xmax()-dfarbox, self.bbox:ymax()-dfarbox,self.bbox:zmax()-dfarbox)
	
	  farbox :pt1(tostring(self.bbox:xmin())..'-dfarbox',tostring(self.bbox:ymin())..'-dfarbox',tostring(self.bbox:zmin())..'-dfarbox')
	  farbox :pt2(tostring(self.bbox:xmax())..'+dfarbox', tostring(self.bbox:ymax())..'+dfarbox',tostring(self.bbox:zmax())..'+dfarbox')
	
	-- farbox :angle_theta(90,360,1)
	farbox :angle_theta('90-elevation_angle','90-elevation_angle',1)
	--farbox :angle_phi(0,360,361)
	farbox :angle_phi(farbox_info.angle.start,farbox_info.angle.stop,farbox_info.angle.steps)
	 farbox :freq(strFreq)
	-- farbox :is_point(true)
	-- farbox :is_outline(true)
	 farbox :is_fill(false)
	farbox :signal(self.psig )	
	 Geometry.AddProbe(farbox)
	----------------------------
end


function RCS:create_grabber(  )
	print( "create grabber")	
	 local cm = Geometry.CGLColorMap:new("Jet0")
	 local  probe = Geometry.TransientFieldsGrabber:new();
	 probe :name("Grabber-xy("..self.name..")");
	 probe :camera(self.camera)
	  probe :save2file(self.save2file);
	 probe :frame_delay(20);
	 probe :size(600,450);
	 
	 local zc = (self.bbox:zmin()+self.bbox:zmax())/2.0
	 probe :add_slice(cm,scalar_e,xy,zc,2,false)
	 -- probe :add_slice(cm,ez,yz,0,2)
	  -- probe :is_fill(false);
	 -- probe:fields_types(ez,hx)
	 probe :stride(20);
	 Geometry.AddTransientFieldsGrabber(probe)
	----------------------------
end



-----------------------------------------------------------------------------
--前处理
function RCS:pre_run ()
end
-----------------------------------------------------------------------------

-----------------------------------------------------------------------------
--后处理
function RCS:post_run ()
	print( "data processing..." )	
	--self.results = RCSPost:new();
	self.results:read_from_model();
	
	print( "get mono RCS..." )	
	self.results:get_mono_RCS();	
	
	print( "plot farfields..." )	
	self:plot_farfield();	
	
	print( "save farfields..." )
	Geometry.find_probe("fbox") :save( self.path .. "/" .. self.name..".txt", "dBsm" )
	
	print( "save geometry..." )
	Geometry.SaveGeometry( self.path .. "/" ..self.name..".gap")
	print( "save summary..." )
	Geometry.save_summary( self.path .. "/" ..self.name.."-calculation-info.txt")
	
	print( "model("..self.name..") run finisted" )
end
-----------------------------------------------------------------------------


-----------------------------------------------------------------------------
--后处理
function RCS:plot_farfield ()
	if self.farfieldView ~= nil then
		self.results.view = self.farfieldView;
		self.results:add_phi_curves(self.results.freq_idx,
							   self.results.theta_idx,
							   self.results.phi_idx,
							   {"RCS-Total"})
	end
end
-----------------------------------------------------------------------------

return RCS;